The Canadian government’s past weather and climate service provides data from weather stations around the whole country. Fortunately, it offers a search by proximity function. Via this service, some sample data of the closest stations to Montreal were downloaded. Each of the data files contains information about the weather stations, including the geographical position (latitude and longitude). These coordinates were plotted on a map (see Figure on the right), and the closest station to the bulk of the data was chosen (station name: MCTAVISH). Only data from this station was used.

Figure: Plot shows all starting stations of a bike trip (red dots), as well as the closest weather stations (blue markers). The closest station in the center is the MCTAVISH weather station (Climate Identifier 7024745, WMO Identifier 71612)

To get a better understanding of the data, the number of hourly bike trips was visualized for the time span between \(2014\) and \(2017\).

The moving average that is shown in the plot (red line) can be interpreted as a baseline model, i.e., the simplest possible model to describe the hourly number of bike rides.

This baseline model explains \(38.8\%\) of the variance \((r^2 = 0.388)\) and has a mean absolute error of \(MAE = 316.2\), which means that on average, the “prediction” for the number of hourly bike rides is wrong by this many bike rides. This includes also winter months with no rides. For a more realistic estimation of model quality, these numbers drop to \(r^2 = .079\) and \(MAE = 510.7\) when only considering the time frame from May to September.

Figure: Number of hourly rides from \(2014\) to \(2017\). Each dot represents the number of trips in one specifc hour. Red line represents a moving average using a window of \(14\) days.

To visualize the relations between the available features, a correlation heatmap is shown on the right. The features are only slightly correlated, with the only exception being temperature and dew point that show an almost perfect (linear) relationship \((r = .93)\).

To avoid problems resulting from this multicollinearity, only temperature was used as a predictor, and dew point was dropped. Despite the fact that that gradient boosting is less influenced by multicollinearity, it might still influence calculations of variable importance [3,4].

Figure: Pearson Correlations between available features in the data.

Table: Model performance measures for different models. \(r^2\) is the amount of variance explained, \(MAE\) stands for mean absolute error. \(train\) and \(test\) specify training and testing set, respectively. For the test set, some performance measures were also re-computed for using only the summer months of the test set (May to September), indicated via the \(summer\) subscripts.

The explained variance of the different models tried ranged from \(r_{test}^2 = 0.860\) to \(r_{test}^2 = 0.941\) for the final model (all in the test set; see table). The hyperparameters for the final model (gradient boosting via scikit-learn) were selected via randomized search using \(4\)-fold cross validation (using only the training set).

The final model fits the data very well, explaining \(94.1\%\) of the variance and exhibiting an average error of \(85.4\) rides over all hourly predictions. This error only increases to \(105.4\) if only summer months are considered.

Figure: The figure shows actual (true) values vs. predicted values for the number of hourly bike trips for the final model. A perfect model would yield predictions that are identical to the true values, i.e., all points would be on the \(45°\) diagonal. This model is relatively close.

The most important features that influenced the prediciton of hourly bike rides were temperature and atmospheric pressure. While the former is easily comprehensible, the latter is best understood as a proxy for precipitation, which was not available from the weather data: low pressure is commonly related to rainy weather [5]. While relative humidity is not as tightly connected to rain [6], it interacts with temperature, e.g., it influences how (high) temperatures are perceived [7].

Further important predictors are the hour of the day and the day of the week, as well as wind direction and speed. How these features influence the predicted number of bike trips is best detailed by specific plots, so-called partial dependence plots (PDP) [1] and individual conditional expectation (ICE) plots [2].

Figure: Most important features in the best-fitting model, gradient boosting via scikit-learn. The plot show the most important features and their relative importance for predicting the hourly number of bike trips. The importance is indicating how well a feature performs in order to reduce the mean squared error for predicting the outcome, averaged over all trees in the ensemble. Values are normed and cannot be interpreted on an absolute scale [8].

Regarding hour of the day, rush hours are clearly associated with higher number of hourly bike rides. There is a peak at at eight o’clock in the morning and a (less pronounced) one in the afternoon. It is also visible that hour of the day interacts with other variables. First, the average PDP line is misleading, as there is a dense area of parallel lines indicating higher bike rides, as well as a high density at zero bike rides (i.e., the average PDP line obscures these two clusters). Second, There are a few lines showing a different pattern, starting to pick up only after the rush hour and steadily rising until the afternoon. This might be the interaction with the day of the week (see section Interactions; see also figure caption for more details about PDP and ICE plots).

Figure: Partial dependence plots (PDP) [[1,2] show the influence of a feature by detailing the effect on the target value (\(y\)-axis) over the whole range of the feature (\(x\)-axis). That is, the thick, dark blue line in the plot shows the main effect of a feature.
In addition, individual conditional expectation (ICE) plots [2] quantify the extent of interactions that exist with other variables. The fine, light blue lines in the plot show predictions at different levels of the feature in question. If they are parallel, there are no interactions with other features. If they are not and show varying slopes compared to the main affect, the feature interacts with other features.
In addition, the plot shows the distribution of the data over the whole feature range as dark blue ticks below the \(y\)-axis. For hour of the day, this is not very interesting (only full hours), but it shows areas with sparse data for other features.

Bike rides also rise with the temperature, the most important feature, once it has reached a certain treshold of about two or three degrees celcius. The rise is relatively steady up to a temperature or about 25 degrees celcius, where the number of hourly rides starts to level off and then to slightly decline. It is also visible (from the spread of the ICE plot lines) that there are massive interactions of temperature with other features.

For relative humidity, the main effect is less pronounced. The general tendency is a slight decrease for higher humidity, with a more pronounced drop around \(95\%\). But the main effect is relatively weak compared to the overall spread of the ICE plot lines, and there are strong interactions with other features for higher relative humidity.

Figure: Partial dependence plots (PDP; thick, dark blue lines) [1] and individual conditional expectation (ICE) plots [2] (thin, light blue lines). For a more detailed explanation, see first occurence of this plot type on the previous page.

The interaction between day of the week and hour of the day is not the strongest interaction in this model (see note on the side), but it serves well to describe the two-feature partial dependence plots (for details, see figure caption). It is clearly visible that during week days (\(0\)-\(4\)), there are peaks in the number of predicted hourly bike rides between \(7\) and \(8\) o’clock in the morning, as well as an even more pronounced peak in the afternoon (higher predicted numbers, as well as a broader peak). The number or predicted bike rides is slightly lower on Monday mornings.

During the weekend (days \(5\) and \(6\)), the distribution is different (i.e., there is an interaction between the two features). There are by far fewer predicted bike rides, also the peak is in the afternoon. Additionally, there are more bike rides during the night, as compared to weekdays.

Figure: Partial dependence plots (PDPs) [1,2] for two features. These plots show the two-way interactions of features (as the number of predicted hourly bike rides, coded as the color in the plot) as a function of the features, in this case day of the week and hour of the day. These plots do not show further, higher-order interactions with other features.

Note: The strength of two-way interactions was determined by estimating a gradient boosting regression model similar to the model described before, but with using main effects and two-way interactions as distinct features.

Figure: Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure … Figure …